/*=================================================================
 *
 * clembek.c	.MEX file corresponding to core routine for ClementsBekkers.M
 * Does fitting of template to a set of data...
 *
 * The calling syntax is:
 *
 *		[pklist, eventlist, crit, scale, cx, nevent] = clembek(template, data, sign, threshold);
 *
 *  You may also want to look at the corresponding M-code, which is commented out
 * in ClementsBekkers.,
 *
 * This is a MEX-file for MATLAB.  
 * Paul B. Manis, Ph.D.
 * 9/13/02 
 *=================================================================*/
/* $Revision: 1.0 $ */
#include <math.h>
#include "mex.h"

/* Input Arguments */

#define	TMPL_IN	prhs[0]
#define	DATA_IN	prhs[1]
#define THRESH_IN	prhs[2]
#define SIGN_IN	prhs[3]

/* Output Arguments */

#define	PKL_OUT	plhs[0]
#define EVL_OUT plhs[1]
#define CRIT_OUT plhs[2]
#define SCALE_OUT plhs[3]
#define CX_OUT plhs[4]
#define N_OUT plhs[5]

#if !defined(MAX)
#define	MAX(A, B)	((A) > (B) ? (A) : (B))
#endif

#if !defined(MIN)
#define	MIN(A, B)	((A) < (B) ? (A) : (B))
#endif

#define PI 3.14159265


//   clembek2( p_pkl, p_evl, p_crit, p_scale, p_cx, p_nout, templ, data, thr, sign, nt, nd); 

static void clembek2(
		   double	*p_pkl, //peaklist (output)
		   double	*p_evl, //eventlist (output)
		   double	*p_crit, // critera at detect (output)
		   double	*p_scale, // scale factor *output)
		   double	*p_cx,		// offset constant (output) 
		   double		*p_nout,	// number of detected events (output)
		   double	*p_templ,	// template (input)
		   double	*p_data,	// data array (input)
		   double	thr,		// threshold for crit detection (input)
		   int	sign,		// sign (pos or neg) for events (input)
		   int		nt,			// number of points in template (input)
		   int		nd			// number of points in the data (input)
		   )
{
	double sume, sume2, sumey, sumy, sumy2, s, c, sse, criteria;
	int i, j, ipk, pkcount;
	double md;


	*p_nout = 0.0;
	pkcount = 0;
	sume = 0.0;
	sume2 = 0.0;
	for (i = 0; i < nt; i++) {
		sume += *(p_templ+i);
		sume2 += (*(p_templ+i) * *(p_templ+i));
	}
	for (i = 0; i < nd-nt+1; i++) {
		if(i == 0) {
			sumy = 0.0;
			sumy2 = 0.0;
			for (j = 0; j < nt; j++) {
				sumy += p_data[j];
				sumy2 += p_data[j] * p_data[j];
			}
		}
		else {
			sumy += (p_data[i+nt-1] - p_data[i-1]);
			sumy2 += ((p_data[i+nt-1] * p_data[i+nt-1])   -  (p_data[i-1]*p_data[i-1]));
		}
		sumey = 0.0;
		for (j = i; j < nt+i; j++) {
			sumey += (p_data[j]*p_templ[j-i]);
		}
		s = (sumey - sume * sumy/nt)/(sume2 - sume * sume/nt);
		c = (sumy - s*sume)/nt;
		/* we don't calculate f as it is not used... */
		sse = sumy2 + (s*s*sume2)+(nt*c*c)-2*(s*sumey + c*sumy - (s*c*sume));
		criteria = s/sqrt(sse/(nt-1));
		*(p_crit+i) = criteria;
		*(p_scale+i) = s;
		*(p_cx+i) = c;
		if(sign == 0 && criteria <= -thr) {
			pkcount += 1;
			*(p_evl+pkcount -1) = i;
			md = 0.0;
			ipk = i;
			for (j=i;j<nt+i; j++) {
				if(p_data[j] < md) {
					md = p_data[j];
					ipk = j;
				}
			}
			*(p_pkl + pkcount - 1) = ipk+1;
		}
		if(sign == 1 && criteria >= thr) {
			pkcount += 1;
			*(p_evl+pkcount-1) = i;
			md = 0.0;
			ipk = i;
			for (j=i;j<nt+i; j++) {
				if(p_data[j] > md) {
					md = p_data[j];
					ipk = j;
				}
			}
			*(p_pkl + pkcount - 1) = ipk+1;
		}

	}

	*p_nout = pkcount; /* now go ahead and update the output count */
}

/* Mex function wrapper for clembek.c */
void mexFunction( int nlhs, mxArray *plhs[], 
		  int nrhs, const mxArray*prhs[] )
     
{ 
    double *p_pkl, *p_evl, *p_crit, *p_scale, *p_cx, *p_nout; 
    double *templ, *data, thr; 
	int nt, nd, sign;
    
    /* Check for proper number of arguments */
    
    if (nrhs != 4) { 
	mexErrMsgTxt("clembek.dll: Four input arguments required."); 
    } else if (nlhs != 6) {
	mexErrMsgTxt("clembek.dll: Requires exactly 6 output arguments."); 
    } 
     
    templ = mxGetPr(TMPL_IN);
	data = mxGetPr(DATA_IN);
	thr = mxGetScalar(THRESH_IN); 
    sign = (int) mxGetScalar(SIGN_IN);
	
/*	 find the number of points in the arrays */
	nt = mxGetN(TMPL_IN);
	nd = mxGetN(DATA_IN);
 
	/* Create matrices for the return arguments */ 
    /* we assume that we cannot detect any more events than points coming in (!) */
	PKL_OUT = mxCreateDoubleMatrix(1, nd, mxREAL); 
	EVL_OUT = mxCreateDoubleMatrix(1, nd, mxREAL); 
	CRIT_OUT = mxCreateDoubleMatrix(1, nd, mxREAL); 
	SCALE_OUT = mxCreateDoubleMatrix(1, nd, mxREAL); 
	CX_OUT = mxCreateDoubleMatrix(1, nd, mxREAL); 
	N_OUT = mxCreateDoubleMatrix(1, 1, mxREAL);
	/* we will return the actual number of valid events in N_OUT */

    /* Assign pointers to the various parameters */ 
	p_pkl = mxGetPr(PKL_OUT);
	p_evl = mxGetPr(EVL_OUT);
    p_crit = mxGetPr(CRIT_OUT);
	p_scale = mxGetPr(SCALE_OUT);
	p_cx = mxGetPr(CX_OUT);
	p_nout = mxGetPr(N_OUT);

   /* Do the actual computations in a subroutine */
    clembek2(p_pkl, p_evl, p_crit, p_scale, p_cx, p_nout, templ, data, thr, sign, nt, nd); 
 // mexPrintf("detected %d events\n", *p_nout);

    return;
    
}


